;;;--------------- solve algbraic eqn with one variable
;;; Problem: solve x in f(x)=0.

;; load streams and range
(library-directories '("." "../../../libs-for-chezscheme"))
(import (range) (stream))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;Bisection method
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Problem: solve x in f(x)=0 in (a b) with f(a)*f(b)<0


;;------ Bisection method
;;; Problem: solve x in f(x)=0, in (a b) with f(a)*f(b)<0.
;;1. i_n => i_{n+1}
(define (next-interval interval f);a interval is a pair (cons x0 x1)
  (let* ([x0 (car interval)]
         [x1 (cdr interval)]
         [x2 (/ (+ x0 x1) 2.0)])
    (if (<= (* (f x0) (f x2)) 0)
        (cons x0 x2)
        (cons x2 x1))))
;;2. interval stream (i0 i1 i2..)
(define (Bisection-stream f init-interval);i0 = init-interval
  (stream-cons init-interval
               (Bisection-stream f
                                 (next-interval init-interval f))))
;;3. solution stream (x0 x1 x2 ..)
(define (Bisection-solutions f a b)
  (stream-map (lambda (interval)
                (/ (+ (car interval) (cdr interval))
                   2.0))
              (Bisection-stream f (cons a b))))

;;;--test: f(x) = x^3 - x - 1, in (1,2)
(let ([f (lambda (x)
           (- (expt x 3) x 1))])
  (for-each (lambda (k)
            (display (stream-ref (Bisection-solutions f 1 2) k))
            (newline))
            (range 10))
  (f (stream-ref (Bisection-solutions f 1 2) 9)))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Newton-Scant method
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;1. x_n => x_{n+1}
(define (next-point/ns x0 x1 f);=>x2
  (- x1 (* (- x1 x0) (f x1)
           (/ 1.0 (- (f x1) (f x0))))))
;;2. solution stream
(define (Newton-Secant-solutions f x0 x1)
  (stream-cons x0
               (Newton-Secant-solutions f x1 (next-point/ns x0 x1 f))))

;;--test: f(x) = x^3 - x -1, in (1,2) has a solution
(let ([f (lambda (x)
           (- (expt x 3) x 1.))])
  (for-each (lambda (k)
            (display (stream-ref (Newton-Secant-solutions f 1 2) k))
            (newline))
            (range 10))
  (f (stream-ref (Newton-Secant-solutions f 1 2) 9)))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Mixed method
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;1. i_n => i_{n+1}
(define (next-interval/mix interval f)
  (let* ([x0 (car interval)]
         [x1 (cdr interval)]
         [x2 (let [(xns (next-point/ns x0 x1 f))
                   (xmid (/ (+ x0 x1) 2.0))]
               (if (and (> xns x0) (< xns x1))
                   xns
                   xmid))])
    (if (<= (* (f x0) (f x2)) 0)
        (cons x0 x2)
        (cons x2 x1))))
;;2. interval stream (i0 i1 i2..)
(define (Mixed-stream f init-interval);i0 = init-interval
  (stream-cons init-interval
               (Mixed-stream f
                             (next-interval/mix init-interval f))))
;; 3. solution stream
(define (Mixed-solutions f a b)
  (stream-map (lambda (interval) (if (< (abs (f (car interval)))
                                        (abs (f (cdr interval))))
                                     (car interval)
                                     (cdr interval)))
              (Mixed-stream f (cons a b))))


;;;----test: f(x) = x^3 - x - 1, in (1,2)
(let ([f (lambda (x)
           (- (expt x 3) x 1))])
  (for-each (lambda (k)
              (display (stream-ref (Mixed-stream f (cons 1 2)) k))
              (newline))
            (range 10))
  (f (stream-ref (Mixed-solutions f 1 2) 9)))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Acceleration method
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Aitken accelerator
(define (Aitken-accelerator strm)
  (let* ([x0 (stream-ref strm 0)]
         [x1 (stream-ref strm 1)]
         [x2 (stream-ref strm 2)]
         [denomi (+ x0 x2 (* -2 x1))])
    (stream-cons (if (= denomi 0.0)
                     x2
                     (- x2 (/ (expt (- x2 x1) 2)
                              denomi)))
                 (Aitken-accelerator (stream-cdr strm)))))
;; stream of streams
(define (stream-of-streams trans s)
  (stream-cons s (stream-of-streams trans (trans s))))
;; super accelerator
(define (Super-accelerator-stream accelerator s)
  (stream-map stream-car (stream-of-streams accelerator s)))
;; The final solutions
(define (Super-Mixed-stream f a b)
  (Super-accelerator-stream Aitken-accelerator (Mixed-solutions f a b)))
(define (Super-Newton-Secant-stream f a b)
  (Super-accelerator-stream Aitken-accelerator (Newton-Secant-solutions f a b)))

;;;----test: f(x) = x^3 - x - 1, in (1,2)
(let ([f (lambda (x) (- (expt x 3) x 1))])
  (for-each (lambda (k)
              (display (stream-ref (Super-Mixed-stream f 1 2) k))
              (newline))
            (range 7)))

;;--- compare four method
(let ([f (lambda (x)
           (- (expt x 3) x 1))])
  (display "Bisection method's error-value: ")
  (display (f (stream-ref (Bisection-solutions f 1 2) 9)))
  (newline)
  (display "Newton-Secant method's error-value: ")
  (display (f (stream-ref (Newton-Secant-solutions f 1 2) 9)))
  (newline)
  (display "Mixed method's error-value: ")
  (display (f (stream-ref (Mixed-solutions f 1 2) 9)))
  (newline)
  (display "Super-Mixed method's error-value: ")
  (display (f (stream-ref (Super-Mixed-stream f 1 2) 9)))
  (newline)
  (display "Super-Newton-Secant method's error-value: ")
  (display (f (stream-ref (Super-Newton-Secant-stream f 1 2) 6)))
  (newline))


